/*
 * HkfModelTest.cpp
 *
 *  Created on: 16 Aug 2011
 *      Author: Allan
 */

#include "HkfModelTest.h"

// C++ includes
#include <iostream>
#include <iomanip>
#include <cmath>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "../../Applications/Equilibrator.h"
#include "../../Utilities/Range.h"

void TestHkfModel_WaterNaCl(double T, double P)
{
	System system({"H2O", "NaCl(aq)"});
	
	Equilibrator equilibrator(system);
	
	cout << setw(15) << left << "mNaCl";
	cout << setw(15) << left << "IonicStrength";
	cout << setw(15) << left << "Act(Water)";
	cout << setw(15) << left << "ActCoeff(NaCl)";
	cout << endl;

	BOOST_FOREACH(double mNaCl, Range(0.1, 6.0, 0.1))
	{
		equilibrator.Reset();
		
		equilibrator.AddCondition("ChargeBalanced");
		equilibrator.AddCondition("FixedMass", "H2O", 55.508);
		equilibrator.AddCondition("FixedTotalMass", "Na[+]", mNaCl);
		equilibrator.AddCondition("FixedTotalMass", "Cl[-]", mNaCl);
		
		SystemState state = equilibrator.Equilibrate(T, P);
		
		const double gNa = state.Get("ActivityCoefficient", "Na[+]");
		const double gCl = state.Get("ActivityCoefficient", "Cl[-]");
		
		const double gNaCl = std::sqrt(gNa * gCl);
		
		cout << setw(15) << left << mNaCl;
		cout << setw(15) << left << state.Get("IonicStrength");
		cout << setw(15) << left << state.Get("Activity", "H2O");
		cout << setw(15) << left << gNaCl;
		cout << endl;
	}
}

void TestHkfModel_WaterCaCl2(double T, double P)
{
	System system({"H2O", "CaCl2(aq)"});
	
	Equilibrator equilibrator(system);
	
	cout << setw(15) << left << "mCaCl2";
	cout << setw(15) << left << "IonicStrength";
	cout << setw(15) << left << "Act(Water)";
	cout << setw(15) << left << "ActCoeff(NaCl)";
	cout << endl;

	BOOST_FOREACH(double mCaCl2, Range(0.1, 3.0, 0.1))
	{
		equilibrator.Reset();
		
		equilibrator.AddCondition("ChargeBalanced");
		equilibrator.AddCondition("FixedMass", "H2O", 55.508);
		equilibrator.AddCondition("FixedTotalMass", "Ca[2+]", mCaCl2);
		equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 2.0*mCaCl2);
		
		SystemState state = equilibrator.Equilibrate(T, P);
		
		const double gCa = state.Get("ActivityCoefficient", "Ca[2+]");
		const double gCl = state.Get("ActivityCoefficient", "Cl[-]");
		
		const double gCaCl2 = std::pow(gCa * gCl * gCl, 1.0/3);
		
		cout << setw(15) << left << mCaCl2;
		cout << setw(15) << left << state.Get("IonicStrength");
		cout << setw(15) << left << state.Get("Activity", "H2O");
		cout << setw(15) << left << gCaCl2;
		cout << endl;
	}
}
